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Introduction 


Recently, significant progress has been made in the handling of singular and nearly- 
singular potential integrals that commonly arise in the Boundary Element Method (BEM) 
[1-3]. To facilitate object-oriented programming and handling of higher order basis 
functions, cancellation techniques are favored over techniques involving singularity 
subtraction. However, gradients of the Newton-type potentials, which produce 
hypersingular kernels, are also frequently required in BEM formulations. As is the case 
with the potentials, treatment of the near-hypersingular integrals has proven more 
challenging than treating the limiting case in which the observation point approaches the 
surface. Historically, numerical evaluation of these near-hypersingularities has often 
involved a two-step procedure: a singularity subtraction to reduce the order of the 
singularity, followed by a boundary contour integral evaluation of the extracted part. 
Since this evaluation necessarily links basis function, Green’s function, and the 
integration domain (element shape), the approach ill fits object-oriented programming 
concepts. Thus, there is a need for cancellation-type techniques for efficient numerical 
evaluation of the gradient of the potential. 

Progress in the development of efficient cancellation-type procedures for the gradient 
potentials was recently presented [2, 4]. To the extent possible, a change of variables is 
chosen such that the Jacobian of the transformation cancels the singularity. However, 
since the gradient kernel involves singularities of different orders, we also require that the 
transformation leaves remaining terms that are analytic. The terms “normal” and 
“tangential” are used herein with reference to the source element. Also, since 
computational formulations often involve the numerical evaluation of both potentials and 
their gradients, it is highly desirable that a single integration procedure efficiently handles 
both. 


Transformation and Method for Near-Hypersingularities 


In order to simplify the presentation and evaluation of the proposed methods, we consider 
the following integral representative of the actual integrals with respect to the 
convergence of numerical integration: 
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where R = |r-r'| is the distance between an observation point at r and source points on 

triangular domains V . Here, the function b is intended to capture the functional 
content of the basis function in scalar form. Following previous approaches [1-5], a 
nearby observation point is projected onto the surface containing the source parent 
triangular element. The parent triangle is then split into three subtriangles about the 
projected observation point. The geometry of a typical subtriangle and its local 
rectangular coordinate system with origin at the projected observation point is shown in 
Fig. 1(a). The general form of the transformed integral over the subtriangle is given by 
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where g(x,y ) is a product of the basis and the Green’s function. In determining a 
suitable transformation, care must be exercised so as not to introduce non-analytic or 
nearly singular functions. A transformation previously investigated [2, 4] can be 
summarized as v = Inf? and u = In [tan (^/2)] . Term s involving odd powers of p , 

which are non-analytic in v, must be handled carefully. In addition, this so-called R 2 
transformation can result in slow convergence when the difference between limits, 
v v -v,, is large. 




Fig. 1 - (a) Subtriangle coordinates and assigned parameters, (b) Subtriangle and disk elements. 

To address these difficulties with the R 2 transformation, we introduce a disk of radius p 

r 0 

centered at the projection point as shown in Fig. 1(b), along with three truncated 
subtriangles. The introduction of this fourth sub-element provides three important 
advantages. First, the azimuthal variation of the gradient can be integrated exactly using 
only two equi-spaced sample points for constant basis functions or with three equi-spaced 
sample points for linear basis functions (completeness order p = 0). The second 
advantage arises from the elimination of the non-analytic terms, p , over the disk. 
Thirdly, as we now show, the remaining radial variation of the integrand can be 
accurately and efficiently integrated. Assuming a basis function of linear order, and 
application of the radial transformation, R = In v , the inner integral over the disk element 
can be written 
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where the functions cir 2 {</>) represent sums of second order terms of the form 

sin 2 (^), cos 2 (f*i), and sin ( </>) cos {(f), and the dependency of R on v has been suppressed for 
notational simplicity. We assume that the source element is sufficiently small that the 
phase term can be expressed accurately as a Taylor series truncated to only a few terms, 
and employ a quadrature rule to numerically integrate the resulting exponential terms. 
We obtain such a rule by solving the system of equations 
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in which J3 = lnji^/z). We can precompute a table of quadrature rules to cover the 
desired range of 8 = pj z . To use the tables, we choose the nearest tabulated entry 
8 = pjz, then modify the disk radius p 0 — > p 0 and corresponding parameter,^ , to 
match the selected table entry. Orders through m ma are integrated exactly over the disk, 
and the insensitivity of the method to the radius, p 0 , is demonstrated in the next section. 

A table accommodating values 0.5 < 8 < 10 5 and disk radius increments of twenty 
percent can be achieved with a table of approximately 70 quadrature rules. We treat the 
remaining truncated subtriangles using the transformations v = Ini? and 

u = In [tan (^/2)] . 


Numerical Results 

The method described above is tested with the source element shown in Fig. 1(b) and for 
observation points at (0, -0.4, z) meters, where z is varied for three cases, 10 ' . 0.01, and 

0. 2 m. The basis function variation is chosen as b - x + y (Fig. 1(b)), and the disk radius 
is set to 0.1 m (except for one case that demonstrates insensitivity to this parameter). For 
reference solutions, x-y integration over the triangle of Fig. 1(b) is achieved with a 
discretization cell size that increases with increasing distance from the projection point; 

1. e., an ^-refinement. For efficient integration of the exponential terms over the disk, we 
choose m rn:a = 2 , which entails a 2-point quadrature rule. Each truncated subtriangle was 

analyzed individually to determine near-optimal ratios for sampling in order to achieve a 
near monotonic convergence for both the normal and tangential gradient components. 
Convergence trends of the normal and tangential components are shown in Figs. 2(a) and 
2(b), respectively, with the number of sample points indicating the combined total for the 
disk and all three subtriangle regions. As seen in Fig. 2(a), accurate integration of the 
normal gradient component is extremely efficient as z tends to zero. The crossed square 
symbols in this figure denote disk quadrature schemes with fewer than six sample points. 
The cross symbols in Fig. 2(a) represent use of a disk radius p 0 = 0.09 , and it can be seen 

that the effects of adjusting the disk radius are minimal. Fig. 2(b) also shows the 
convergence results for the potential using the same quadrature rules developed for the 
gradient. Clearly, the same rule can be used for both. 



Conclusion 


A technique for efficient numerical evaluation of hypersingular kernels is demonstrated. 

It is also shown that weaker singularities may be handled efficiently by the same 
technique. As the projection of the observation point approaches an edge of the source 
element, the efficiency of the technique is reduced. Treatment of this case remains as 
future work. Additional work also includes extending the technique to handle projections 
of the observation point falling outside the source element. 



Number of sample points Number of sample points 

Fig. 2 - Convergence of components (a) normal gradient and (b) tangential gradient and potential. 
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